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ABSTRACT 



Shipboard fires have been the bane of mariners since man's earliest attempts to sail 
the sea. Understanding the behavior of fire in an enclosed space, such as those found 
on today's modern seagoing vessels, will greatly enhance the mariner's ability to combat 
or prevent them. In a joint effort between the Naval Postgraduate School and the 
University of Notre Dame a computer code has been developed to model a full scale fire 
in a closed compartment. The code uses a finite volume formulation to obtain numerical 
solutions to the unsteatdy, three dimensional conservation equations of mass, momen- 
tum and energy. Included are the effects of turbulence, strong buoyancy, surface radi- 
ation and wall conduction. The code gives velocities, pressure, temperatures, and 
densities throughout the field. 

This thesis applies that computer code to the U.S. Navy's full scale fire test chamber 
at Naval Air Warfare Center, China fake, California. Advanced computer graphics 
techniques, including color contouring and three dimensional vector field plotting have 
been applied to make output data more informative. It is hoped that someday this 
model could provide a useful tool for naval architects in the design of a fire safe ship, 
and a cost effective means for development 'evaluation of new firefighting equipment and 
techniques. 



c.i 

THESIS DISCLAIMER 

T he reader is cautioned that computer programs developed in this research may not 
have been exercised for all eases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 



A. BACKGROUND 

Fire aboard ship is one of the professional mariner's worst fears. With no place to 
go to escape the heat and smoke the mariner must fight the blaze or face the loss of his 
vessel. Even though L .S. merchant marine and naval personnel receive some of the best 
fire fighting training available, annual losses to shipboard fires can easily run into the 
hundred million dollar range. Ship down time, equipment repair replacement, personnel 
injuries and casualties all contribute to these costs and result in the degradation of our 
merchant marine and naval forces. To minimize these losses it is imperative that the 
phenomena of fire be studied and understood especially in the relatively closed environ- 
ment found aboard ship. It is only through study and understanding that adequate 
means may be developed to prevent or mitigate the devastating effects of a shipboard 
fire. 

Fire is a complex phenomena whose study is equally complex, requiring the com- 
bined knowledge of a variety of fields including fluid dynamics, heat and mass transfer, 
and combustion. Research into the mechanics of fire and development of methods to 
predict it's behavior will aid engineers in reducing the probability of its ignition and 
minimizing its effects once ignited. [Ref. 1] 

It is the complexity of fire which makes its study so difficult, especially in a 
shipboard environment. Small scale fire studies have proven inadequate for predicting 
the behavior of large scale fires. Therefore, full scale studies are generally the only 
means for conducting realistic experimentation. However, such full scale experimenta- 
tion can be very costly and dangerous. Shipboard fires often occur in fully enclosed, 
airtight spaces where pressures can build up during the fire. These spaces may have re- 
stricted accesses, contain electronic equipment, fiammables and or toxics, all of which 
add to the expense and danger of full scale experimentation. Additional costs will also 
be incurred as air quality and emission standards are stiffened across the country. 

To study the phenomena of shipboard fire, the U.S. Navy has built a test facility at 
the Naval Research Laboratory in Washington, DC known as Fire-1. This facility is 
essentially a large cylindrical pressure vessel with spherical endcaps intended for use in 
full scale modeling of fires inside submarines and'or in closed airtight compartments 
aboard ships. Missile attacks against British warships during the Falklands crisis and the 
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Iraqi missile attack on the USS STARK in the Persian Gulf in 1987 prompted the L'.S. 
Navy to build a second fire test facility at the Naval Air Warfare Center (N’AWC) in 
China Lake, CA. The purpose of this facility is to study the clTects of fire in a vented 
compartment |Ref. 2, pp. 7-8]. It is this facility which is modeled in this thesis. 

A less expensive, and less dangerous alternative to full scale fire experimentation is 
the use of computer modeling techniques, 'fhe development of high speed, high capacity 
computers has allowed researchers to model thoroughly the complex fire phenomena 
and predict fire behavior without the expense of full scale testing. A properly developed 
computer model, validated against actual full scale test data, can provide a less expensive 
and safe alternative to full scale experimentation. Furthermore, the inherent flexibility 
of computer modeling allows it to be used on increasingly more complex geometries. 
Someday this may lead to modeling of entire ships during the design phase to identify 
areas particularly susceptible to fire, or for the accurate prediction of the effectiveness 
of new firefighting techniques. 

B. COMPUTER MODELING 

Field modeling uses finite difference techniques to subdivide the volume being 
studied, in this case the simulated shipboard compartment, into small, finite volume el- 
ements. Using initial conditions specified by the user and the finite difference forms of 
the equations for conservation of mass, momentum, and energy, values of temperature, 
pressure, velocity and density are calculated for each of the individual volumetric ele- 
ments at discrete time intervals. Additional modeling of physical effects such as radi- 
ation, turbulence and wall conduction are included to increase the validity of the 
simulation. The enormous number of calculations necessary for this type of modeling 
mandates the use of large amounts of computer memory and high speed processors. 

The basis for this thesis is provided by a large number of previous research projects. 
One of the earliest successful models of this type was developed by Aziz and Heliums 
|Ref. 3] at Rice University in Houston, Texas. They expressed the Navier-Stokes 
equations in terms of vorticity and vector potential, then solved the resulting three di- 
mensional finite difference equations using a combination of the alternating direction 
method and successive over-relaxation (SOR). Later work by Mallinson and de Vahl 
Davis | Ref. 4], Morrison and Tran |Ref. 5]; Chan and Banerjee [Ref. 6]; Ozeo, Fujii, Lior 
and Churchill [Ref. 7|; and a host of others have all expanded on the. use of finite dif- 
ference techniques to model convective heat transfer in closed compartments. 
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In the late 1970's and early 1980's, R.G. Rehm and ll.R. Baum began developing 
equations to describe the buoyant flow induced by large scale fires (Refs. 8, 9, 10 and 
1 1] Work done at the University of Notre Dame used a two dimensional finite difference 
field model to predict velocities, temperatures and smoke concentrations in aircraft cabin 
fires (Ref. 12 and 13]. The development of a two dimensional model of transient, natural 
convection cooling was developed, and experimentally verified, by Xicolette el al (Ref. 
14] using a semi-implicit upwind differencing scheme and global pressure corrections. 
Still more studies (Refs. 15, 16, 17, 18 , & 19] have utilized finite difference methods to 
solve non-linear, three dimensional partial differential equations for rectangular enclo- 
sures. 

C. THE FIRE TEST FACILITY 

A full scale test facility has been constructed at the Naval Air Warfare Center 
(NAWC) at China Lake, CA. This facility, designed to simulate full scale shipboard 
compartments, consists of three chambers (see Figure 1 on page 4 and Figure 2 on page 
5]. The main chamber measures 20 feet by 20 feet by 10 feet and is vented to the at- 
mosphere through an opening in a side wall, in a manner intended to simulate the pen- 
etration of a missile or other projectile. The second chamber measures 15 feet by 15 feet 
by 10 feet is located adjacent to the main chamber to the east, while the third chamber, 
also measuring 15 feet by 15 feet by 10 feet, is located on top of the main chamber. 
These secondary chambers are intended to provide a means to study the vertical and 
horizontal heat transmission rates. (Ref. 2, pp. 7-10] 

The entire test structure is constructed of 3/8 inch thick steel bulkheads 'walls rein- 
forced by 5 inch 1 -beams on 5 foot centers, and 1/2 inch thick steel decks reinforced 
underneath by 12 inch 1-beams on 5 foot centers. Additional support was provided for 
the overhead in the main chamber by several 5 inch 1-beam columns. T hese columns 
were intended to eliminate any sagging and/or distortion of the overhead as it is repeat- 
edly cycled from ambient conditions to the extreme temperatures expected during the 
missile fuel tests. Access to all three compartments is through hatches which open to 
the outside of the structure and are kept closed during all testing. There are no openings 
between any two adjacent chambers. The ventilation opening is located on the north 
face of the main chamber to simulate an impact hole, such as a missile strike (see 
Figure 3 on page 6 and Figure 4 on page 7], (Ref. 2, pp. 10-12] 

Instrumentation is provided within the test chamber to measure radiation heat flux, 
total heat flux, compartment bulk pressure, wall temperatures and compartment gas 
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Figure 1. Fire test chamber at NAWC, China Lake, CA. 
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Figure 2. Plan view of fire test chamber at NAWC, China Lake, CA. 
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Figure 3. North/south elevations of test chamber at NAWC, China Lake, CA. 
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Figure 4. East/west elevations of test chamber at NAWC, China Lake, CA. 
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temperatures. Thermocouples were arranged in vertical arrays of 10 thermocouples. 
The lowest thermocouple in each array was placed 6 inches above the deck with subse- 
quent thermocouples spaced at 1 foot intervals ending 6 inches below the ceiling [see 
Figure 5 on page 9). [Ref. 2, pp. 14-15) 

D. THE COMPUTER PROGRAM 

This computer model is a joint project undertaken by the Naval Postgraduate 
School and the University of Notre Dame. It represents a low cost, safe alternative to 
full scale fire testing. With proper modifications, and properly validated by full scale 
experimentation, this program should provide a valuable tool for testing the effectiveness 
of fire mitigation techniques and evaluation of new ship designs. 

Work on this program began at Naval Postgraduate School in 1986 when Nies [Ref. 
20] used a cartesian coordinate system to model the cylindrical spherical geometry of the 
FI RE-1 test facility. In 19S7 Raycraft [Ref. 2 1 ] modified the program to use a 
cylindrical spherical coordinate system. She also expanded the scope of this project by 
writing a companion program to calculate view factors and account for surface radiation 
effects. Using both of Raycraft's programs, Houck [Ref. 22] further expanded the scope 
of the project to account for the internal ventilation capabilities of the FI RE- 1 facility. 
Most recently, in 1991, McCarthy [Ref. 23 1 used advanced computer graphics techniques 
to provide a more accurate pictorial representation of the unique geometry of the 
FIRE-1 facility and demonstrated the advantages of using full color displays to represent 
the three dimensional isotherm and velocity vector field profiles. This development 
greatly enhanced the presentation of the program output. McCarthy was also the first 
researcher at NPS to incorporate Raycraft's radiation/view factor program as subpro- 
grams of the main computer code. 

This thesis returns to the cartesian coordinate system used by Nies [Ref. 20], to 
model the Navy's newest fire test facility at the Naval Air Warfare Center (NAWC) at 
China Lake, CA. It also advances McCarthy's [Ref. 23] work by applying the enhanced 
graphics capabilities of the NCAR Graphics software, developed by the National Center 
for Atmospheric Research [Refs. 24 and 25], for output presentation. 
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Figure 5. Instrumentation of fire test chamber at NAWC, China Lake, CA. 
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II. THE NUMERICAL MODEL 



A. THE GOVERNING EQUATIONS. 

The model used in this code is built on the fundamental laws of conservation of mass 
(continuity), energy and momentum. One example of these equations is presented by 
Patankar [Ref. 26, pp. 11-17). As in that work, this model neglects the pressure work 
and viscous dissipation terms of the energy equation due to the low flow velocities ex- 
pected. The volumetric heat generation term will be used to account for heat input by 
the fire. Patankar's version of the momentum equations [Ref. 26, p. 14, eq. 2.11] in- 
cluded both a body force and an extra viscous force. This model assumes a newtonian 
fluid therefore the extra viscous force is neglected and the only body force present is due 
to gravitational effects acting in the negative z-direction. Thus the body force in the 
z-momentum equation is set equal to — pg where g is the local gravitational constant; 
the negative sign indicates that g is acting in the negative z-direction; and p is the fluid 
density. 

Xies [Ref. 20, pp. 16-38] expanded the equations presented by Patankar into the 
three dimensional, cartesian coordinate system used in this model. Following the de- 
velopment pattern of Doria |Ref. 18, pp. 4-7) and making the assumption that air is a 
homogeneous gas of constant composition (reactive variations due to the fire are neg- 
lected), Nies [Ref. 20, pp. 16-17] incorporated the equations of state for an ideal gas with 
constant specific heats in the form: 



P= pRT 



and 



h = c p (T-T re/ ) 

where P, T and p represent the bulk pressure, temperature and density inside the control 
volume, R is the gas constant for air, h is the specific enthalpy, c p is the specific heat and 
T re/ is a suitably chosen reference temperature. 

Again, following the procedures developed by Patankar [Ref. 26] and expanded by 
Doria |Ref. 18], Nies [Ref. 20 , pp. 21-26] goes on to place these six equations in their 
non-dimensional, integral forms. After subdividing the test chamber into a large, but 
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finite number of control volumes, finite difference techniques are used to solve for the 
six unknowns of temperature, pressure, density and the three components of velocity. 

B. THE CONTROL VOLUME. 

Each control volume, or cell, is defined by the nodal point contained at its center. 
Values for temperature, density and pressure are calculated at this central point (desig- 
nated P) and are then assumed to hold for the entire cell. A secondary, or staggered grid 
system, is used to determine velocities. This staggered grid is offset from the main grid 
by one half the length of the cell [sec Figure 6 on page I2|. As explained in McCarthy 
[Ref. 23, pp. 18-19) and Patankar [Ref. 26, pp. I 15-120], use of this staggered grid alle- 
viates two problems: first, since velocities are dependent upon pressure differentials, the 
staggered grid allows pressures to be determined at a more frequent interval thereby re- 
ducing the error associated with larger separations between nodal points; second, the 
stagger increases stability by decreasing and or eliminating unrealistic, oscillatory veloc- 
ity fields when adjacent velocities arc used to satisfy the continuity equation. 

Since Patankar's method uses primitive variables, in lieu of the stream function and 
vorticity, special attention must be paid to the coupling of equations through the pres- 
sure terms. An iterative process is used to calculate the pressure in each cell, then a local 
pressure correction is calculated to ensure that continuity is satisfied. Additional dis- 
cussions of this correction are included in both Patankar [Ref. 26, pp. 120- 128| and Doria 
[Ref. 18, pp. 26-32]. A global pressure correction was described by Nicolette, el al [Ref. 

1 4 1 for addressing net energy changes within the closed compartment. This correction 
has also been incorporated into this model. 

As stated in McCarthy [Ref. 23, p. 19], forcing convergence on a non-linear set of 
equations, like those describing fluid flow, can be difficult. A variety of schemes have 
been developed, each with its own set of strengths and weaknesses. I his model applies 
iterative techniques to a solution scheme known as "QUICK", or Quadratic Upstream 
Interpolation for Convective Kinematics, developed by Leonard [Ref. 19). With the ac- 
curacy of a central differencing scheme and the stability of the convective diffusion terms 
in upwind differencing, the QUICK scheme estimates values and gradients of the trans- 
port variables at the cell faces, flic utility of this method was demonstrated by fl.Q. 
Yang [Ref. 27] when he used QUICK to solve the coupled momentum and energy 
equations for three dimensional flow in a tilted rectangular enclosure. 

In addition to the center cell described above, neighboring cells will be used in vari- 
ous computations. To keep track of these cells a standardized nomenclature must be 
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Figure 6. One Dimensional Basic and Staggered Computational Cells. 



adopted. Assuming the central cell, cell P, to be located at the point (i.j,k) where i.j and 
k correspond to the standard coordinated axis x. y and z respectively, then each neighbor 
may be described as follows (NOTE: the following directions are for standardization of 
nomenclature only and do not necessarily correspond to compass directions as shown 
in Figure 1 on page 4): 

East (i -r 1 ,j,k) 

West (i-l.j.k) 

North (i.j a- l,k) 

South (i,j-l,k) 

Front (i.j,k+- 1) 

Back (i.j,k-l) 

The nodal point in each direction is designated by the capital letter corresponding 
to the direction (i.e. E, W, N, S. F, B), while the boundary between cell P and its 
neighbors is designated by a lower case letter corresponding to that direction (i.e. the 
boundary between cell P and cell N is designated n, and the boundary' between P and F 
is 0 - A typical cell array on a rectangular grid is shown in figure Figure 7 on page 13. 




As discussed previously, each P node is used to determine values of density, pressure 
and temperature which are them applied to its entire cell. Velocities are determined at 
the cell faces using the staggered grid arrangement described above. 

C. DISCRETIZATION OF THE CONSERVATION EQUATIONS. 

With that brief background, the integral forms of the conservation equations devel- 
oped Nies [Ref. 20, pp. 25-26] may now be discretized. As discussed by Nies [Ref. 20, 
pp. 26-38] maximum stability and accuracy of the model is achieved by using three dif- 
ferent finite differencing schemes. Forward differencing is used for timewise 
discretization, central differencing is used to discretize the diffusion terms, and the 
QUICK scheme is used to discretize the convective terms. Use of these techniques to 
discretize the governing equations was discussed in detail by Nies [Ref. 20, pp 26- 38] 
and will not be repeated here. The following finite difference forms of the governing 
equations are taken from his work. 

1. The Continuity Equation. 

n AjrArA? 

(p P - Pp) — + (G e - GJAyAz + (G n - G s )AzAjc + (G f - G b )AxAy = S mp 
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where p P is the density at the current time step; p° p is the density at node P at the previous 
time step; At is the incremental time step; Ajc, Ay, and A 2 are the cell dimensions in the 
indicated direction; and S p m is the residual mass term. The mass flux (G) terms are de- 
fined as: 



n ( Pe + Pp) 

t' M e 



G w = 



( Pp + P ir) 



Gn = 



(Px + Pp) 



^ ( Pp + Ps ) 

G ‘ = 2 r ' 



^ (Pp + Pp) 

G f = j w f 



n (Pp + Pb ) 

Gfr 2 

It should be noted that in all of these equations p refers to the density; u, v, and 
vv are the three velocity components in the x, y, and z-directions respectively; uppercase 
subscripts denote values at the indicated nodal point; and lowercase subscripts denote 
values at cell faces. Also, in order for continuity to be satisfied in this closed system, the 
residual mass term ( S mp ) would equal zero. However, due to the approximation inherent 
in the numerical scheme, we will be satisfied if S mp tends toward zero as determined by 
comparison to an arbitrarily small threshold value. 

2. The Energy Equation. 



IL 



o A-vAyAz 
+ Pp' 



At 



h P = ll A h E + H Aw h\v+ UaJis + H A s h S + H. A r h F + h aJ^b + H s P 



where h is the specific enthalpy at the current time step; h° is the specific enthalpy at the 
previous time step; and the coefficients (H) are defined as: 



H a e = 



(I G e \-G e ) 



AyAz + 



AyAz 
Re [ Pr l J e Ax 
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,, (I G w \ + G w ) 
u a w = j A - v ’ Az + 



(\G H \-G H ) t t 
U An = j AzA - v + 



„ (\ G s\+Gs) A A 

II As = j AzA - v + 



(Kyi - Gy) 

II a — A* Ay 4- 



„ (1^1+^) AA 
II Ab = Ax Ay + 



] \ AyAz 


RC'Pr, 


) w Ax 


1 


\ AzAx 


Re l Pr l 


F 4 y 




\ AzAx 


Re,Pr, 


) s 4 y 




) A-vAy 


RC'Pr, 


y a z 




\ AxAy 


ReJX 


N 

<1 

*C> 



"a, = " A E +"a w + II A , + II As + II A f + II A B 
II s p = Pph Q pAxAyAz 

Also used above arc the turbulent Reynolds number (Re r ) and the turbulent Prandtl 
number ( Pr ,). These arc defined as 



Re t = 



Po^pII 

*V/ 



Pe/fpO 

b 

k eff 

where 11 is the characteristic length (defined as the height of the test chamber in our 
model); p and c p represent the density and specific heat of air, while the subscript 0 in- 
dicates that the properties are to be evaluated at the initial conditions which existed prior 
to ignition of the fire; Uq is the reference velocity (set equal to 1.0 ft/sec); and n rff and 
A e// are effective values of viscosity and conductivity as defined by Nies |Ref. 20, p. 39-40). 
3. The Momentum Equations. 

As stated by Nies [Ref. 20, p. 31] the momentum equations are more complex 
than the previous equations because of the use of the staggered grid and the addition 
of the shear stress terms. Staggered grids are determined by shifting the main grid one 
half cell in the negative direction along one axis at a time. Maintaining the nomencla- 
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turc established for the basic, centered cells, the central node of a cell staggered in the 
x-dircction becomes w, the central node in a cell staggered in the y-dircction becomes s, 
and the central node in a cell staggered in the z-direction becomes b. Likewise, the faces 
of the staggered cell are designated by the capital letters representing the basic nodes 
through which they pass. Thus, in the staggered grid P always represents the positive 
face along the axis where the shift (stagger) is being evaluated. Thus, the x-momentum 
equation becomes: 



0 AwAyAz 

"b Pw 



u w ~ ^e tl e + d WM U ww + 



.YvvTVvv ~b ^Sw^Sw ~b s^Fw^Fw ~b ^Bw^Bw ~b 



where 






(I G P | - G,) , ( Re, V 1, . 
^ 1 : AvAz 



Ax 



A ww 



(I C W ,J +G h ,,) ( Re, U 1 4 

H : AvAz 



Ax 



^ Xw 



( D., 



(\G nw \-G nH ) | v Re t 



Ay 



AzAx 






( | G sw | + G sw ) + Re r 



( D ^ 



Ay 



AzAx 



A Fw 



(\G/ W \ - G fw ) y Re t 



( Rp 



+ ■ 



Az 



Ax Ay 



^Bw 



(I G bw \+G bw ) ( Re, )bw 1 A 
■ H : Ax Ay 



Az 



^ e ~b -^ww "b ^ AW ~b "b ~b ^Bw 



16 



S„ = p!«I 



AjrAt’Az 

A/ 



-(Pp- P w)A>Az 



+ ^»H’) 



1 

Ke, 



AyAz 

Ajc 



- (u lv - u w J 



1 

Ke, 



Vi' Vi' 



AyAz 

A.v 



+ 



( V AV 



* A’hu) 



1 






Az 



+ 



0‘V - 






Ay 



and 




u P = 



u u = 



“> = 






0 . 0 
Pp + P» 



“e + “u 
2 

»u- + «»v 
2 

»AV + »H- 
2 

»5 m' + 



«Fw + 
2 

d” 

2 



G/> — PpW/> 
G\y = p nil w 



G 



m v 



f Pa + Pp , Pa'm 7 + P ir 

[_ 2 V " + 2 lW J 



2 



17 



